A MULTIRESOLUTION SPACE-TIME ADAPTIVE SCHEME FOR THE BIDOMAIN 

MODEL IN ELECTROCARDIOLOGY 

MOSTAFA BENDAHMANE^, RAIMUND BURGER^, AND RICARDO RUIZ BAIER^ 

Abstract. This work deals with the numerical solution of the monodomain and bidomain models of elec- 
trical activity of myocardial tissue. The bidomain model is a system consisting of a possibly degenerate 
parabolic PDE coupled with an elliptic PDE for the transmembrane and extracellular potentials, respec- 
tively. This system of two scalar PDEs is supplemented by a time-dependent ODE modeling the evolution 
of the so-called gating variable. In the simpler sub-case of the monodomain model, the elliptic PDE re- 
duces to an algebraic equation. Two simple models for the membrane and ionic currents are considered, 
^"^ the Mitchell-Schaeffer model and the simpler FitzHugh-Nagumo model. Since typical solutions of the bido- 

^^ main and monodomain models exhibit wavefronts with steep gradients, we propose a finite volume scheme 

enriched by a fully adaptive multiresolution method, whose basic purpose is to concentrate computational 
vN effort on zones of strong variation of the solution. Time adapt ivity is achieved by two alternative devices, 

namely locally varying time stepping and a Runge-Kutta-Fehlberg-type adaptive time integration. A series 

^^ of numerical examples demonstrates that these methods are efficient and sufficiently accurate to simulate 

"^Vj the electrical activity in myocardial tissue with affordable effort. In addition, an optimal threshold for 

\^ discarding non-significant information in the multiresolution representation of the solution is derived, and 

• the numerical efficiency and accuracy of the method is measured in terms of CPU time speed-up, memory 

^iQ compression, and errors in different norms. 
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1. Introduction 



^ 1.1. Scope. The obvious difficulty of performing direct measurements in electrocardiology has motivated 

C^ wide interest in the numerical simulation of cardiac models. In 1952, Hodgkin and Huxley [ ] introduced 

^V the ffist mathematical model of wave propagation in squid nerve, which was modified later on to describe 

several phenomena in biology. This led to the first physiological model of cardiac tissue [33] and many others. 
Among these models, the bidomain models firstly introduced by Tung [ ], is one of the most accurate and 
complete models for the theoretical and numerical study of the electric activity in cardiac tissue. The 
bidomain equations result from the principle of conservation of current between the intra- and extracellular 
^^ domains, followed by a homogenization process (see e.g. [ , 14, 28]) derived from a scaled version of a 

cellular model on a periodic structure of cardiac tissue. Mathematically, the bidomain model is a coupled 
system consisting of a scalar, possibly degenerate parabolic PDE coupled with a scalar elliptic PDE for the 
transmembrane potential and the extracellular potential, respectively. These equations are supplemented 



H by a time-dependent ODE for the so-called gating variable, which is defined at every point of the spatial 

computational domain. Here, the term "bidomain" refiects that in general, the intra- and extracellular 
tissues have different longitudinal and transversal (with respect to the fiber) conductivities; if these are 
equal, then the model is termed monodomain models and the elliptic PDE reduces to an algebraic equation. 
The degenerate structure of the mathematical formulation of the bidomain model is essentially due to the 
differences between the intra- and extracellular anisotropy of the cardiac tissue [ , ] . 

The bidomain model represents a computational challenge since the width of an excitation front is roughly 
two orders of magnitude smaller than the long axis of a human-size right ventricle. This local feature, along 
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2 BENDAHMANE, BURGER, AND RUIZ 

with strongly varying time scales in the reaction terms, produces solutions with sharp propagating wave fronts 
in the potential field, which almost precludes simulations with uniform grids. Clearly, cardiac simulations 
should be based on space- (and also time-) adaptive methods. 

It is the purpose of this paper to develop a fully adaptive multiresolution (MR) scheme with locally 
varying space-time stepping (LTS) and adaptive time step control by means of a Runge-Kutta-Fehlberg 
(RKF) method. These strategies are of different nature, but do not exclude each other; rather, they may 
be combined to obtain a potentially more powerful method (as is suggested e.g. in [ ]; however, herein 
we do not pursue that approach). We furthermore address the deduction of an optimal threshold value for 
discarding non-significant data, which permits to achieve significant data compression. Previous experience 
with degenerate parabolic equations and reaction-diffusion systems [ , 7, 8, 9, 36, 37] suggests that the 
MR device should provide an efficient tool for solving the bidomain equations, and in this same spirit, we 
construct the corresponding extension of the MR method with the novel application to the bidomain and 
monodomain models in mind. 

The efficiency of the MR method is a consequence of the fact that at each time step, the solution is encoded 
with respect to a MR basis corresponding to a hierarchy of nested grids. The size of the details determines 
the level of refinement needed to obtain an accurate local representation of the solution. Therefore, an 
adaptive mesh is evolved in time by refining and coarsening in a suitable way, by means of a strategy based 
on the prediction of the displacement and creation of singularities in the solution. 

We apply the MR approach to an explicit finite volume (FV) method in each time step. Since the 
computational effort required for integrating a system of equations for one time step is usually substantially 
higher for an implicit scheme when compared to explicit schemes, implicit schemes may be less efficient than 
explicit ones, especially when the overall number of time steps is large (see e.g. [ ]). 

1.2. Related work. To further put this work into the proper perspective, we first mention that stan- 
dard theory for coupled parabolic-elliptic systems (see e.g. [IC]) does not apply naturally to the bidomain 
equations, since the anisotropics of the intra- and extracellular media differ and the resulting system is of 
degenerate parabolic type. Colli Franzone and Savare [ ] present a weak formulation for the bidomain 
model and show that it has a structure suitable for applying the theory of evolution variational inequalities 
in Hilbert spaces. Bendahmane and Karlsen [ ] prove existence and uniqueness for the bidomain equations 
using the Faedo-Galerkin method and compactness theory for the existence part, and Bourgault et al. [ ] 
prove existence and uniqueness for the bidomain equations by first reformulating the problem as a single 
parabolic PDF, and then applying a semigroup approach. 

From a computational point of view, substantial contributions have been made in adaptivity for cardiac 
models. However, the approach presented herein differs to the best of our knowledge from other adaptive 
approaches in the literature. These alternative techniques include adaptive mesh refinement (AMR) (e.g., 
[11, 4i]), adaptive finite element methods using a posteriori error techniques (see, e.g., [ ]) or multigrid 
methods applied to finite elements. Furthermore, Quan et al. [35] present a domain decomposition approach 
using an alternating direction implicit (ADI) method. With respect to time adaptivity, Skouibine et al. [39] 
present a predictor-corrector time stepping strategy to accelerate a given finite differences scheme for the 
bidomain equations using active membrane kinetics (Luo-Rudy phase II). Cherry et al. [ ] use local time 
stepping, similar to the method introduced in the germinal work of Berger and Oliger [ ], to accelerate a 
reference scheme. Parallelized versions of part of the methods mentioned above are presented, for example, 
by Colli Franzone and Pavarino [15] and Saleheen and Ng [ ]. 

MR schemes for hyperbolic partial differential equations were first proposed by Harten [ ]. We refer to 
the work of Miiller [ ] for a survey on MR methods, see also Chiavassa et al. [i.]. As stated above, the 
idea behind the MR method is to accelerate a reference discretization scheme while controlling the error. 
In the context of fully adaptive MR methods [ ], the mathematical analysis is complete only in the case 
of a scalar conservation law, but in practice, these techniques have been used by several groups (see e.g. 
[2, 21, 30, 31, 37]) to successfully solve a wide class of problems, including applications to multidimensional 
systems. For more details on the framework of classical MR methods for hyperbolic partial differential 
equations, we also refer to Cohen et al. [ ] and Dahmen et al. [ ]. 
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1.3. Outline of the paper. The remainder of this paper is organized as fohows. In Section 2, the bidomain 
and monodomain models of cardiac tissue are introduced. The general bidomain model can be expressed as 
a coupled system of a parabolic PDE and an elliptic PDE plus an ODE for the evolution of the local gating 
variable, while the monodomain model, which arises as a particular sub-case of the bidomain model, is defined 
by a reaction-diffusion equation, which is again supplemented with an ODE for the gating variable. Section 3 
deals with the construction of an appropriate FV method for the solution of both the parabolic-elliptic system 
and the reaction-diffusion equation arising from the bidomain and monodomain models, respectively. Next, 
in Section 4 we develop the MR analysis used to endow the reference finite FV schemes with space adaptivity. 
More precisely, in Section 4.1 we introduce the wavelet basis underlying the multiresolution representation 
with the pertinent projection operator. In Section 4.2, the prediction operator and the detail coefficients are 
introduced. Small detail coefficients on fine levels of resolution may be discarded (this operation is called 
thresholding), which allows for substantial data compression. In Section 4.3, we recall the graded tree data 
structure used for storage of the numerical solution, and which is introduced for ease of navigation. In 
Section 4.4 we outline an error analysis, similar to that conducted in [ , , ] and motivated by the rigorous 
analysis of Cohen et al. [ ], which justifies the choice of a reference tolerance £r. In turn, this quantity 
determines the comparison values si used for the thresholding operation at each level / of multiresolution. 
The basic goal is to choose the threshold values in such a way that the resulting multiresolution scheme has 
the same order of accuracy as the usual finite volume scheme. 

In Section 5 we address two strategies for the adaptive evolution in time of the space-adaptive MR scheme, 
namely the locally varying time stepping (LTS, Section 5.1) and a variant of the well-known Runge-Kutta- 
Fehlberg (RKF, Section 5.2) method. Finally, in Section 6 we present numerical examples putting into 
evidence the efficiency of the underlying methods. Some conclusions that can be drawn from the paper 
about the effectiveness of our methods and statement of possible further extensions to our research are given 
in Section 7, and in the Appendix we present a brief description of the LTS and general MR algorithms. 

2. The macroscopic bidomain and monodomain models 

The spatial domain for our models is a bounded open subset ^^ C M^ with a piecewise smooth boundary 
dCt. This represents a two-dimensional slice of the cardiac muscle regarded as two interpenetrating and 
superimposed (anisotropic) continuous media, namely the intracellular (i) and extracellular (e) tissues. These 
tissues occupy the same two-dimensional area, and are separated from each other (and connected at each 
point) by the cardiac cellular membrane. The quantities of interest are intracellular and extracellular electric 
potentials, u^ = Ui{x^ t) and Uq = Uq{x^ t), at (x, t) G VLt := ^ x (0, T). Their difference v = v{x^ t) := u^ — Uq 
is known as the transmembrane potential. The conductivity of the tissue is represented by scaled tensors 
Mi(x) and Me(x) given by 

M,(x) = cT*I + (cTj.-(T*)a,(x)a;r(x), 

where aj = crj(x) G C"^(M^) and a* = crj(x) G C"^(R^), j G {e,i}, are the intra- and extracellular conductivi- 
ties along and transversal to the direction of the fiber (parallel to a^(x)), respectively. 

For fibers aligned with the axis, Mi(x) and Me(x) are diagonal matrices: Mi(x) = diag(<j|, a*) and 
Me(x) = diag( (7^,(7*). When the so-called anisotropy ratios cr|/crf and cl/al are equal, we are in the case of 
equal anisotropy, but generally the conductivities in the longitudinal direction 1 are higher than those across 
the fiber (direction t); such a case is called strong anisotropy of electrical conductivity. When the fibers 
rotate from bottom to top, this type of anisotropy is often referred to as rotational anisotropy. 

The bidomain model is given by the following coupled reaction-diffusion system [^zu, ^o]: 

PCmdtV - V • (Mi(x)Viii) + Plioniv, w) = 0, 
pCmdtV + V • (Me(x)Viie) + PIion{v, w) = /app, (2.1) 

dfW — H{v, w) = 0, (x, t) G ^T- 

Here, Cm > is the so-called surface capacitance of the membrane, /3 is the surface-to-volume ratio, and 
w{x,t) is the so-called gating or recovery variable, which also takes into account the concentration variables. 
The stimulation currents applied to the extracellular space are represented by the function /app = /app(x,t). 
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The functions H{v^ w) and /ion ('^7 u)) correspond to the fairly simple Mitchell- Schaeffer membrane model [29] 
for the membrane and ionic currents: 

RmCmVooiv/Vp) Rm \Vpr]2 V^r]i J 

where the dimensionless functions r]oo{s) and Woo{s) are given by r]oo{s) = ^3 + (^4 — ^3)^(5 — 7^5) and 
'^00 (-5) = H{s — 7^5), where H denotes the Heaviside function, R^ is the surface resistivity of the membrane, 
and Vp and 771, . . . , 775 are given parameters. A simpler choice for the membrane kinetics is the widely known 
FitzHugh-Nagumo model ['^^, '^'^], which is often used to avoid computational difficulties arising from a large 
number of coupling variables. This model is specified by 

H{v,w) = av — bw^ Iion{v^w) = —X(w — v{l — v){v — 6))^ (2-3) 

where a, 6, A and are given parameters. 

We rewrite (2.1) equivalently in terms of v and Uq as the strongly coupled parabolic-elliptic PDE-ODE 
system (see e.g. ['^, ^^]): 

Pc^dtV + V • {Me{x)\/Ue) + PIion{v. w) = /app, (2.4a) 

V • ((Mi(x) + Me(x))V^Xe) + V • (Mi(x)V^) = /app, (2.4b) 

dtw- H{v,w) = 0, {x,t)eVtT- (2.4c) 

We utilize zero fiux boundary conditions, representing an isolated piece of cardiac tissue: 

[mj{x)Vuj) • n = on St := 5(^ X (0,T), j G {e,i}, (2.5) 

and impose initial conditions (which are degenerate for the transmembrane potential v): 

v({)^x) = vq{x)^ w({)^x) = wq{x)^ x^Vt. (2.6) 

For the solution v of the bidomain model, we require the initial datum vq to be compatible with (2.5). 
Therefore, if we fix both 1x^(0, x), j G {e, i} as initial data, the problem may become unsolvable, since the 
time derivative involves only v = u^ — Uq (this is also referred as degeneracy in time). Thus, we impose the 
compatibility condition 

/ u^{x, t)dx = for a.e. t e (0, T). (2.7) 

Jn 

In the case that Mi = AMe for some constant A G M, the system (2.1) is equivalent to a scalar parabolic 

equation for the transmembrane potential v^ coupled to an ODE for the gating variable w. This parabolic 

equation is obtained by multiplying the first equation in (2.1) by 1/(1 + A), the second by A/(l + A) and 

adding the resulting equations. The final monodomain model can be stated as follows: 



fSCmdtV - V • ( :r^^^ ) + /^^ion(^, W) = :^^^/app, 



Mi ^ \ ,, , , A 

^^V.j+/?/ion(.,-) = ^^^.pp, ^^g^ 

dfW — H{v^ w) = 0, (x, t) e fir- 



This model is, of course, significantly less involved and requires substantially less computational effort than 
the bidomain model, and even though the assumption of equal anisotropy ratios is very strong and generally 
unrealistic, the monodomain model is adequate for a qualitative investigation of repolarization sequences 
and the distribution of patterns of durations of the action potential [ ] . 

We assume that the functions Mj, j G {e,i}, /ion, and H are sufficiently smooth so that the following 
definitions of weak solutions make sense. Furthermore, we assume that /app G L'^{Qt) and Mj G L'^{Q) 
and Mj^ • ^ ^ C'mICP foi" a.e. x e ft, for all ^ G M^, j G {e,i}, and a constant Cm > 0. For later reference, 
we now state the definitions of a weak solution for the bidomain and the monodomain model, respectively. 

Definition 2.1. A triple u = (v^Uq^w) of functions is a weak solution of the bidomain model (2.4)-(2.6) if 
v^Uq G L'^{0,T;H^{Q)), w G C([0,T], L^(1]))^ (2.7) is satisfied, and the following identities hold for all test 
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functions ^^^l^^i G P([0,T) x Ct): 

Pcm / vo{x)(p{0,x)dx^ // \(3cmvdt(p-Me{x)Vue-V(p^(3Iion^\dxdt= // h^^Lpdxdt, 

JQ J JQt J JQt 

// |-(Mi(x) +Me(x))Viie • V?A-Mi(x)Vv- VV^jdxdt = // h^^ifdxdt, 

— / wo{x)^{0^x) dx — // wdtidxdt = II H^dxdt. 

Definition 2.2. A pair u = (v^w) of functions is a weak solution of the monodomain model (2.8) if 
V G L'^{O^T] H^{Q)), w G C([0,T], L^(r^))^ and the following identities hold for all test functions (f^£, ^ 
D([0,T) xCl): 

Pcm / vo{x)(p{0,x)dx 

JQ 
+ / / \ PCmVdt^p + /3/ion^ - :— — -Mi VV • V(/P ^ dx dt = :— — - / / /appV^ ^^ dt, 



n^l 1 + A ^J 1 + A^^^^ 

— / wo{x)^{0^x) dx — // wdtidxdt = II H^dxdt. 

JQ J JQt J JQt 

3. The reference finite volume scheme 

To define the FV scheme for approximating solutions to the bidomain equations (2.4), we follow the 
framework of [ ] . An admissible mesh for Q is formed by a family T of control volumes (open and convex 
polygons) of maximum diameter h. (From the next section on, we will use only Cartesian meshes, however 
the following description is in a more general setting.) For all K G T^ xk denotes the center of K^ N{K) is 
the set of neighbors of K^ Eint{K) is the set of edges of K in the interior of T, and ^ext(^) the set of edges 
of K on the boundary dft. For all L G N{K), d{K^ L) denotes the distance between xk and xl, cfk.l is the 
interface between K and L, and rjK^L {vk.cti respectively) is the unit normal vector to cf^^l {o' G £e:^t{K)^ 
respectively) oriented from K to L (from K to dft^ respectively). For oil K G T^ \K\ stands for the measure 
of the cell K. The admissibility of T implies that Q = UxeTK, Kr\L = 0ifK^LGT and K ^ L^ and 
there exist a finite sequence {xK)KeT for which xJ^xZ is orthogonal to crK,L' 

Now, consider K G T and L G N{K) with common vertices {a£^K,L)i<£<i with / G N\{0} and let Tk,l 
(T^^ for a G Se^t{K)^ respectively) be the open and convex polygon with vertices (x^, ^l) {xr^ respectively) 
and {ai^K,L)i<i<i' Notice that Q can be decomposed into Q = lJKer{{^LeN{K)TK,L) U {^aeSe^tiK)TK,a))' 

For all K G T^ the approximation VhUh of Vu is defined by 

[ul - uk)Vk,l if ^ G Tk,l^ 



\/hUh{x) := { \Tk,l 

To discretize (2.4)-(2.6), we choose an admissible discretization of Qt, consisting of an admissible mesh 
of Vt and a time step size At > 0. For example, we could choose A/" > as the smallest integer such that 
N/Si, > T, and set f^ := nAt for n G {0, . . . , N}. However, for reasons stated further below, we assume that 
the time step size At is determined anew for each iteration, and present the scheme as it applies to advance 
the solution from T to r+^ := t^ + At. 

On each cell K e T, (positive definite) conductivity tensors are defined by 



M,,. 



j,K = Tj^i Mj{x)dx, j G {e,i}. 
Let Fj^K,L be an approximation of 

Mj{x)\/Uj • r]K,Ldj 



I 
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for j e {e, i}, and foiKeT and L G N{K), let 



M, 



J,K,L 



-^ / Mj{x)dxr]K,L 



\Mj^kVk,l\, j e {e,i}. 



Here |-| stands for the Euclidean norm. The diffusive fluxes 'M.j{x)\/uj • r]K,L on (Jk,l are approximated by 



'^ 0-K,L 



U — X Jk^ U ' — U ' T<r 

^ \crK,L\^Uj{ycr) • (Mj^kVk.l) = WK,L\Mj^K,L^Uj{y^) • —^ ^ \(Tk,l\Mj^k,l-^ ^, 

where y^ is the center of (Jk,l and Uj^o- is an approximation oiujijjcr)^ j G {e, i}. The resulting approximation 
of fluxes is consistent [ ]. In addition, the scheme should be conservative. This property enables us to 
determine the additional unknowns Uj^^^ and to compute the numerical fluxes on internal edges: 

Fi,K,L = dlK,L^0-^{ui,L - u^,k) ULe N{K), (3.1) 



where we define 



^* .^ Mj^k,lMj^l,k ..^ J-. 



while we discretize the zero-flux boundary condition by setting 

Fj,K,a = (3.2) 

on boundary edges. Here, d{K^ ^k,l) and (i(L, (Jk,l) are the distances from xk and xl to (Jk,l^ respectively. 
We define cell averages of the unknowns H{v^ w) and /ion('^, u))'- 



.n+l 



^^^' •" Aij]^ /„ J^H{v{x,t),w{x,t))dxdt, III]k-=^^J^^ I^Iion{v{x,t),w{x,t))dxdt, 
and of the given function /app- 



1 r^' f 

The computation starts from the initial cell averages 

VK = n^ vo{x)dx, '^K = n^ wo{x)dx. (3.3) 

1^1 Jk 1^1 Jk 

We now describe the finite volume scheme employed to advance the numerical solution from t^ to t^+^, 
which is based on a simple explicit Euler time discretization. Assuming that at t = t^, the quantities u'^ x^ 
j G {e, i}, 1'^ = {u^^x ~ '^e k)^ ^^^ '^K ^^^ known for all i^ G T, we compute the values of these cell averages 
u]^x\ j e {e, i}, vi^^ = (2i[;+^ - <+^) and <+^ at t = T+i from 

Pc^m"^^^^^ E <K,LJ^«L-<K)+/5|^|/ron,K = l^l4%,K, (3.4) 

LeN{K) ^ ' ^ 

E j^ { Kk,l + dl,K,L) Kr - <%') + dlK,L K+' - «r ') } = 1^14-pp.i., (3.5) 

LeN{K) ^ ' ^ 

ii^ r^^'^;^^ -wgK=o- (3.6) 

The order in which these equations are used to advance is explicitly stated in Algorithm 1 in Section 6. 

The boundary condition (2.5) is taken into account by imposing zero fluxes on the external edges as in 
(3.2): 

4if..J^^«,L-«",K)=0 iovaG£,^t{K), je{e,i}, n = 0,1,2,.... (3.7) 
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and the compatibility condition (2.7) is discretized via 

^ \K\uIk = {), n = 0,l,2, 



KeT 



Analogously, a FV method for the monodomain model (2.8) is given by determining vectors v^ and 
w'^^ foiKeT and n = 0, 1, 2, . . . such that for all i^ G T, we start from the initial data given by (3.3), 
and use the following formulas to advance the solution over one time step: 

^.n+l _ n 1 l/T I A 



LeN{K) 






A FV method for a slightly different version of the bidomain equations is analyzed in [']. In that paper, 
the authors prove existence and uniqueness of solutions to an implicit FV scheme, and provide convergence 
results. On the other hand, following [ , 18, 23], we prove in [ ] existence and uniqueness of approximate 
solutions (that is, well-definedness of the scheme) for an implicit version of the scheme considered herein, 
and show that it converges to a weak solution in the sense of Definition 2.1, under a mild condition that 
limits the time step size At (a CFL-type condition is, however, not imposed in [ ]). Moreover, as in [ ], 
we may deduce that in the case of Cartesian meshes, the explicit version of the FV method utilized herein, 
(3.3)-(3.7) , is stable under the CFL condition 

At < /i(2max(|/r„„,^| + |4"pp,^|) + Ah-^ max(|Mi,^| + iMe.^l))"'- (3.8) 

Notice that the values of /j^n,K ^^^ ^2pp,K depend on time. However, while /app is a given control function 
for our model and therefore max^er l^rpp k\ ^^^ assumed to be bounded, the quantity I^^^ ^ is not bounded 
a priori for arbitrarily large times. Consequently, in our computations, we evaluate the right-hand side of 
(3.8) after each iteration at t = t^, and use (3.8) to define the time step size At to advance the solution from 

r to r+^ = r + At. 

4. MULTIRESOLUTION AND WAVELETS 

4.1. Wavelet basis. Consider a rectangle which after a change of variables can be regarded as Q = [0, 1]^. 
We determine a nested mesh hierarchy Aq C • • • C A/,, using an uniform dyadic partition of Q. Here each 
grid A^ := {V(ij),z}(i,j), with (i, j) to be defined, is formed by the control volumes at each level V(i,j),^ '-= 
2~^[i,i + 1] X [j,jf + 1], z, j G Ii = {0, . . . , 2^ — 1}, / = 0, . . . ,L. Here, / = corresponds to the coarsest 
and I = L to the finest level. The nestedness of the grid hierarchy is made precise by the refinement sets 
'^{i,j),i = {2(^,j) + e}, e G £^ := {0,1}^, which satisfy #M(^ij)^i = 4. For x = (xi,X2) G V^ij^^i the scale 
box function is defined as 

^(i,j),K^) •= pfT \Xv^i,j)A^) = 2^^X[o,i]2(2^xi -i,2^X2 -j), 

and the averages of any function u{-^t) G L^{Q) for the cell V(i,j),z may be expressed equivalently as the 
inner product U(^ij^^i := ('W, ^(i,j),z)Li(n)- We are now ready to define the following two-scale relation for cell 
averages and box functions: 

El^r,Z+l| ~ _ ^ Y^ ~ 

IT/, , |V^r,^+l - 7 2^ V^(2i+p,2j+g),^+l, 



KiJ),l 






yHl)l\ 4 
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which defines a projection operator, which ahows us to move from finer to coarser levels. For x G 'i^2(i,j)+a,^+i 
with a G £^, we define the wavelet function depending on the box functions on a finer level 

7 Y^O-2/ i\a-e~ V^ ^ r=2(i,j)+a,^ + l | ^ -. \a-e ~ 



s.eE rGA4^- .^ - ^^ii^3)M 






The number of related wavelets is #Mi^ij)^i — 1 = 3. Since r • e G {0,1,2} for r, e G £^, we have for 
instance that 

'^(i,i),(l,0),^ = ^ (^2(i,j),^+l + ^2(i,j) + (0,l),^+l - ^2(i,j) + (l,0),^+l - ^2(i,j) + (l,l),^+l) • 

Doing this for all e G £^* := £^ \ {(0, 0)} yields an inverse two-scale relation (see [ ]), namely 

eeE 

This equation is related to the concept of stable completions [ ]. Roughly speaking, the L^-counterparts of 
the wavelet functions {ilj(^ij^^i}ij^i^ form a completion of the L^-counterpart of the basis system {(^(^ .,) J^^^^/^ , 
and this determines the existence of a biorthogonal system. 

4.2. Detail coefficients. For e G £^*, we introduce the details^ which will be crucial to detect zones with 
steep gradients: d(^ij^^^^i := {u^ipi^ij^^^^i). These detail coefficients also satisfy a two-scale relation, namely 

dii,j),.,i = 4 Yl (-l)""^2(^,,)+a,^+i. (4.2) 

2(i,j)+aGAl(i,,-),z 

An appealing feature is that we can determine a transformation between the cell averages on level L and the 
cell averages on level plus a series of details. This can be achieved by applying recursively the two-scale 
relations (4.1) and (4.2); but we also require this transformation to be reversible: 

^iiJ),W= Yl 3[iJ),r^r,h S^J) •= {^([i/2]+ri,[i/2]+r2),J^^,^,e{-5,...,0,...,5}' ('^•^) 

where Sj. -^ is the stencil of interpolation or coarsening set, g\. .^ ^ are coefficients, and the tilde over u in 
the left-hand side of (4.3) means that this corresponds to a predicted value. 

Relation (4.3) defines the so-called prediction operator, which allows us to move from coarser to finer 
resolution levels. In contrast to the projection, the prediction operator is not unique, but we will impose 
two constraints: to be consistent with the projection, in the sense that the prediction operator is the right 
inverse of the projection operator, and to be local, in the sense that the predicted value depends only on 
Sl- -y For sake of notation, in our case we may write (4.3) as 

where ei, 62 G {0, 1} and 

s s 

Qx '= 2^ In {Ui^i^rij),l - U(i-n,j),l) i Qy '= Z^ ^P {^ihj+p),l ~ ^iiJ-p),l) ' 
n=l p=l 

s s 

Qxy '= /_^ln/_^%{^ii+n,j+p),l — ^{i+nj-p),l — '^{i-nj+p),l ^ ^{i-nj-p),l) • 
n=l p=l 

Here the corresponding coefficients are 71 = — ^ and 72 = j|g (see [ ]). 

From [19] we know that details are related to the regularity of a given function: if u is sufficiently smooth, 
then its detail coefficients decrease when going from coarser to finer levels: 

where r = 2s + 1 is the number of vanishing moments of the wavelets. This means that the more regular u is 
over V(ij),^, the smaller is the corresponding detail coefficient. In view of this property, it is natural to attain 
data compression by discarding the information corresponding to small details. This is called thresholding. 



MULTIRESOLUTION SCHEME FOR THE BIDOMAIN MODEL 9 

Basically, we discard all the elements corresponding to details that are smaller in absolute value than a 
level-dependent tolerance £/, 

Given a reference tolerance £r, which is determined by means of an error analysis (see Section 4.4), we 
determine ei by 

ei = 22^^-^)£R. (4.4) 

4.3. Graded tree data structure. We organize the cell averages and corresponding details at different 
levels in a dynamic graded tree: whenever an element is included in the tree, all other elements corresponding 
to the same spatial region in coarser resolution levels are also included, and neighboring cells will differ by at 
most one refinement level. This choice guarantees the stability of the multiscale operations [ ]. We denote 
by root the basis and by node an element of the tree. In two space dimensions, a parent node has four sons, 
and the sons of the same parent are called brothers. A node without sons is called a leaf. A given node 
has s' = 2 nearest neighbors in each spatial direction, called nearest cousins^ needed for the computation 
of the fluxes of leaves; if these nearest cousins do not exist, we create them as virtual leaves. The leaves 
of the tree are the control volumes forming the adaptive mesh. We denote by A the set of all nodes of the 
tree and by C{A) the restriction of A to the leaves. We apply this MR representation to the spatial part of 
the function u = {v^ Uq^w)^ which corresponds to the numerical solution of the underlying problem for each 
time step, so we need to update the tree structure for the proper representation of the solution during the 
evolution. To this end, we apply a thresholding strategy, but always keep the graded tree structure of the 
data. Once the thresholding is performed, we add to the tree a safety zone, so the new tree may contain the 
adaptive mesh for the next time step. The safety zone is generated by adding one finer level to the tree in 
all possible positions without violating the graded tree data structure. This device, first proposed by Harten 
[26], ensures that the graded tree adequately represents the solution in the next time step. Its effectiveness 
depends strongly on the assumption of finite propagation speed of the singularities. 

Note that the fluxes are only computed at level / + 1 and we set the ingoing flux on the leaf at level / 
equal to the sum of the outgoing fluxes on the leaves of level / + 1 sharing the same edge 

Fii+l,j),l^ii,j),l = ^(2i+l,2j),^+l^(2i+2,2j),^+l + ^(2i+l,2j+l),^+l^(2i+2,2j+l),/+l- (4.5) 

It is known that this choice decreases the number of costly flux evaluations without loosing the conservativity 
in the flux computation, and this represents a real advantage when using a graded tree structure, see e.g. 
[6(] for more details. This advantage is lost for a non-graded tree structure, for which fluxes for leaves on 
an immediately finer level are not always available. 

The data compression rate [8, 9] 77 := A/'/(2~*^^+-^W + #>C(A)) is used to measure the improvement in 
data compression. Here, A/" is the number of elements in the full finest grid at level L, and #>C(A) is the 
size of the set of leaves. We also measure the speed-up V between the CPU time of the numerical solution 
obtained by the FV method and the CPU time of the numerical solution obtained by the MR method: 
V := CPUtimepv/CPUtimeMR. 

4.4. Error analysis of the multiresolution scheme. Using the main properties of the reference FV 
scheme, such as the CFL stability condition and order of approximation in space, we can derive the optimal 
choice for the threshold parameter sr for the adaptive MR scheme. We can decompose the global error 
between the cell average values of the exact solution vector at the level L, denoted by u^^ = ('U^, u^^^^, '^i^), 
and those of the MR computation with a maximal level L, denoted by u^j^, into two errors 

ll^ex ^MRII -^ ll^ex ^FV|| ^ ||^FV ^MR||- 

The first error on the right-hand side is called discretization error and the second perturbation error. Using 
this, the CFL condition (3.8), and the fact that h = \ft\^^'^2~^; we obtain that if the so-called reference 
tolerance for the numerical computations in Section 6 is set to 

2-(a+2)L 
£R = C y r y ^, (4.6) 

\n\ max(|/ion,K| + 2|/app,K|) + |^p/'22+^ max(|M,,K| + {M^^kI 
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then we can expect the discretization error and the perturbation error to be kept at the same order (see 
[8, 13, 37] for more details). 

To measure errors between a reference solution iiex and an approximate solution umr^ we will use L^- 
errors: Cp = Wu"^^ - i^mrIIp' P = 1^ 2, oo, where 

V ^ ' iid,i)eCiA) 
Here u^^- ■ ^ stands for the projection of the reference solution onto the leaf (i, j, /). 




5. Time-step accelerating methods 



5.1. Local time stepping. We employ a version of the locally varying time stepping strategy introduced 
by Miiller and Stiriba [31], and summarize here its principles. The basic idea is to enforce a local CFL 
condition by using the same CFL number for all scales, and evolving all leaves on level / using the local time 
step size 

Ati = 2^-^ At, / = L- 1,...,0, 

where At = AtL corresponds to the time step size on the finest level L. This strategy allows to increase 
the time step for the major part of the adaptive mesh without violating the CFL stability condition. The 
synchronization of the time stepping for the portions of the solution lying on different resolution levels will 
be automatically achieved after 2^ time steps using Ati. To additionally save computational effort, the the 
tree is updated only each odd intermediate time step 1, 3, . . . , 2^ — 1, and furthermore, the projection and 
prediction operators are performed only on scales occupied by the leaves of the current tree. For the rest 
of the intermediate time steps, we use the current (old) tree structure. For the sake of synchronization and 
conservativity of the fiux computation, for coarse levels (scales without leaves), we use the same diffusive 
fiuxes and sources computed in the previous intermediate time step, because the cell average on these scales 
are the same that in the previous intermediate time step. Only for scales containing leaves, we compute 
fiuxes in the following way: if there is a leaf at the corresponding edge and at the same resolution level /, 
we simply perform a fiux computation using the brother leaves, and the virtual leaves at the same level if 
necessary; and if there is a leaf at the corresponding cell edge but on a finer resolution level / + 1 {interface 
edge), the fiux will be determined as in (4.5), i.e., we compute the fiuxes at a level / + 1 on the same edge, 
and we set the ingoing fiux on the corresponding edge at level / equal to the sum of the outgoing fiuxes on 
the sons cells of level / + 1 (for the same edge). To always have at hand the computed fiuxes as in (4.5), we 
need to perform the locally varying time stepping recursively from fine to coarse levels. 

5.2. A Runge-Kutta-Fehlberg method. In order to upgrade the FV scheme described in Section 3 to 
at least second order so that the second-order spatial accuracy is effective, we utilize an RKF method [ ], 
which, apart from providing the necessary accuracy, also allows an adaptive control of the time step. For 
our models, we consider a vector-valued RKF method, i.e., u = {v,Uq,w) and its time-discretized form at 
step 771, denoted by u^. For ease of discussion, we assume that the problem is written as 9^u = A{t, u). 

We use two Runge-Kutta methods, of orders p = 3 and p — 1 = 2 



where 



^m+l 


= u"^ + biK.1 + 62/^2 + ^3^3, 


^m+l 


= u'^ 


+ biKi + b2K2 


: + ^3/^3, 




^1 : 


= AU(t^,u^), 














^2 : 


= AtAit"^ + C2 At, 


U^ + 


a2i^i). 










^3 : 


= AU(t^ + C3 At, 


U- + 


asi^i + a32^2), 






]s corresponding 


to the RK3(2) method 


are C2 


= CL21 = 1, 


C3 = 


1 

2 



(5.1) 



<^31 — <^32 — 4, 

5-L = 62 = |, 63 = |, bi =1)2 = ^^ and 63 = 0. These values yield an optimal pair of embedded TVD-RK 
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Time 


V 


V 


L-^— error 


L^— error 


L^ -error 


t = 0.0 ms 




170.22 


4.31 X 10-"^ 


2.47 X 10-"^ 


3.99 X 10-"^ 


t = 1.5 ms 


27.81 


37.56 


4.97 X 10-"^ 


1.96 X 10-"^ 


4.63 X 10-"^ 


t = 3.5 ms 


26.47 


29.89 


5.23 X 10-"^ 


4.05 X 10-"^ 


4.82 X 10-"^ 


t = 4.5 ms 


31.41 


28.12 


7.48 X 10-"^ 


4.29 X 10-"^ 


5.31 X 10-"^ 


t = 5.5 ms 


30.62 


24.70 


1.04 X 10-^ 


6.20 X 10-"^ 


6.79 X 10-"^ 



Table 1. Example 1 (Monodomain model): Corresponding simulated time, CPU ratio V, 
compression rate rj and normalized errors for v^ using a MR method. 



Time [ms] V 77 Potential L^— error L^— error L^— error 

t = 0.1 13.74 19.39 V 3.68 x lO""^ 8.79 x 10"^ 6.51 x 10""^ 

Ue 2.01 X 10-"^ 6.54 X 10-^ 5.22 x lO""^ 

t = 0.5 21.40 17.63 v 4.06 x lO""^ 9.26 x 10"^ 6.83 x lO""^ 

2.79 X 10-"^ 8.72 x 10"^ 5.49 x lO""^ 

A ^^ -X r^—A -t ^r- -, r^—A ^ r^ r^ -x r^ — A 



Ue. 



t = 2.0 25.23 17.74 v 4.37 x lO""^ 1.25 x lO""^ 6.88 x 10""^ 

u^ 3.48 X 10-"^ 9.44 x 10"^ 6.11 x 10""^ 

t = 5.0 26.09 16.35 v 5.29 x lO""^ 1.94 x lO""^ 7.20 x 10""^ 

Ue 4.15 X 10"^ 1.06 X 10"^ 6.32 x 10"^ 

Table 2. Example 2 (bidomain model, one stimulus): Corresponding simulated time, CPU 
ratio V, compression rate rj and normalized errors. 



methods of orders two and three. The truncation error between the two approximations for u^+^ is estimated 

by 

P 

4id := u"^+' - U-+' = Y.(k - h)R, ^ {Atr, ^oid := ||4id||oo. (5.2) 

i=l 

Then we can adjust the step size to achieve a prescribed accuracy (^desired in time. The new time step is 
determined by Atnew = ^^oid I ^desired /^oidT^^ with p = 3. To avoid excessively large time steps, we define 
a limiter function S{t) := {Sq — 5min) exp(— t/At) + 5min, where we choose Sq = 0.1 and Smin = 0.01. The 
new time step At^ew is then defined as 

^^ ^ fAtold|(^desired/(^old|'/^ if |(Atnew - Atold)/Atold | < ^5(t, Atold), .^ 3. 

I ^S{t, Atoid) Atoid + Atoid otherwise. 

Notice that Atnew is the time step size for computing u^+^. More details on the RKF scheme and its 
implementation can be found in [8, 20]. 

6. Numerical Examples 

We will present three test cases showing the efficiency of the previously described methods in capturing the 
dynamical evolution of electro-physiological waves for both the monodomain and bidomain models. Since we 
are dealing with multicomponent solutions, we emphasize that a single mesh is used to represent the vector 
of relevant variables. In the bidomain model, the anisotropics, mesh structures, and the size of the problem 
cause the sparse linear system corresponding to (3.5) to be ill-conditioned. This system needs to be solved in 
each time step, which is done by the Cholesky method. Before presenting the numerical results, we provide 
further details on the implementation of the numerical schemes. 
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Figure 1. Example 1 (monodomain model): Numerical solution for v, measured in [mV] 
(left) and leaves of the corresponding tree at times (from top to bottom) t = 1.5 ms, t = 
3.5 ms. 



6.1. Implementation Issues. The following algorithm shows how the solution u"^+^ = {v^Uq^w)^~^^ is 
obtained in each time step is obtained for each time step 

Algorithm 1 (General method). 

(1) Assume that u^ , u^ , v^ and w^ are known (at time t^). 

(2) Solve the ODE 

dtw — H{v^ w) = 0^ X G Q^ 

approximately for f^ < t < t'^'^^ with initial condition w'^ and data v'^ , i.e., compute w"^^^ using 
(3.6). 

(3) Solve the parabolic PDE 

pC^dfV + V • (Me(x)Viie) + Phon{v, w) = /app, X G 1^, 

(Me(x)V^ie) •n = on SI] 
approximately for t^ < t < t^~^^ , with v{t^) = v^ and w{t^) = w^, i.e., calculate v'^^^ using (3.4). 
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Figure 2. Example 1 (monodomain model): Numerical solution for v, measured in [mV] 
(left) and leaves of the corresponding tree at times (from top to bottom) t = 4.5 ms and 
t = 5.5 ms. 



(4) Solve the elliptic problem 

V- ((Mi(x)+Me(x))Viie)+V- (Mi(x)V^) 



^app? 



X eft 



{Mj{x)\/Uj) • n = on dft, j G {e, i} 

approximately for f^ < t < t'^'^^ with v{t'^) = v^ and UQ{t'^) = u^, i.e., determine u^^^ by solving 
the linear system (3.5). 

This algorithm structure is usually preferred for systems involving parabolic and elliptic equation, since 
it explicitly isolates the solution of the elliptic problem from the rest of the computations [ ]. 

For multicomponent solutions, there are many possible definitions for a scalar detail <^(i,j),/ that is cal- 
culated from the details of the components (see a brief discussion in [ ]). To guarantee that the refine- 
ment and coarsening procedures are always on the safe side, in the sense that we always prefer to keep 
a position with a detail triple containing at least one component above the threshold (4.4), we will use 
dY- -N , = minify. .^ ,, d^^ .^ .,df- -^ ,| and d^- -^ , = maxjdy. .^ ,, d^^ .^ j,df. -^ ,| for the refinement and coarsen- 
ing procedures, respectively. In practice, the details introduced in Section 4.2 are computed simply as the 
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Figure 3. Example 2 (bidomain model, one stimulus): Initial condition for the extracellular 
potential Uq^ and leaves of the corresponding tree data structure. 



differences between the "exact" and the predicted value: 






In (4.6) a stands for the accuracy order of the FV reference scheme, and numerical experience gives 
a = 1.09. To choose an acceptable value for the factor C, a series of computations (not completely shown 
here) with different tolerances are needed in each case, prior to final computations. We basically choose the 
largest available candidate value for C such that the same order of accuracy as that of the reference FV 
scheme (same slope) is maintained and the data compression rate rj and the speed-up V are maximized (see 
Figure 6). In [ ] we give a detailed description of the multiresolution algorithm and the LTS algorithm for 
systems of reaction-diffusion equations. For sake of completeness, we provide versions of these algorithms 
adapted to the bidomain model of electrocardiology in the Appendix. 

6.2. Example 1. For this example, we consider the simple monodomain model (2.8) with homogeneous 
Neumann boundary conditions. The ionic current and membrane model is determined by the FitzHugh- 
Nagumo membrane kinetics (2.3), with a = 0.16875, b = 1.0, A = —100 and = 0.25. The computational 
domain is the square 1] = [0, 1 cm]^, and the remaining parameters are c^ = l.OmF/cm and P = 1.0 cm~^. 
The units for v^w are mV. We consider in (2.8) (1 + A)~^Mi := diag(7, 7) with 7 = 0.01. The respective 
initial data for v and w are 



v^{x,y) 



1 



1 



mV, Wo = mV. 



1 + exp(-50(x2 + 7/2)1/2 _ 0.1) 
After 4 ms, an instantaneous stimulus is applied in (xo,^o) = (0.5 cm, 0.5 cm) to the membrane potential v 



A 



1 + A 



^app 



ImV if {x - xo)^ + (^ - yo)'^ < 0.04 cm^, 
mV otherwise. 



In this example, we use L = 10 resolution levels, Af = 262144 elements in the finest level, a tolerance of 
£r = 1 X 10~^, and we compute normalized errors by comparison with a reference solution obtained with 
a fine mesh calculation with JV = 1024^ = 1048576 control volumes. The time evolution is made using a 
first-order explicit Euler scheme. Plots of the numerical solution with the corresponding adapt ively refined 
meshes at different times are shown in Figures 1 and 2. 
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Figure 4. Example 2 (bidomain model, one stimulus): Numerical solution for transmem- 
brane potential v and extracellular potential Uq in [mV], and leaves of the corresponding 
tree data structure at times t = 0.1 ms and t = 0.5 ms. 



As can be seen from Table 1, the normalized errors are controlled to be of the same order of the reference 
tolerance sr. We also see that the MR algorithm is efficient: we have high rates of memory compression and 
speed-up. 
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Figure 5. Example 2 (bidomain model, one stimulus): Numerical solution for transmem- 
brane potential v and extracellular potential Uq in [mV], and leaves of the corresponding 
tree data structure at times t = 2.0 ms and t = 3.5 ms. 



6.3. Example 2. In Examples 2 and 3, we present computational results for the simulation of the bidomain 
equations. We consider a computational domain Q = [0,5cm]^, and the parameters in (2.2) and (2.4) (after 
[22, 39, 43, 44]) are given by the membrane capacitance Cm = 1.0 mF/cm , the intracellular conductivity in the 
principal axis a\ = 6Q~^cm.~^, the remaining intracellular conductivity a* = 0.6(7~^cm~^ (corresponding 
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Figure 6. Example 2 (bidomain model): data compression rate rj (left), speed-up factor V 
(middle) and L^-errors for different scales L and values of sr (right). The simulated time 
is t = 2.0 ms. 



to an anisotropy ratio of 10), the extracellular conductivities al = 241] ^cm ^ and a* = 12 ^^ ^cm ^ 
(corresponding to an anisotropy ratio of 2), the surface-to- volume ratio /3 = 2000 cm~^, the surface resistivity 
Rm = 2x lO'^l^cm^, Vp = 100 mV, r]i = 0.005, r]2 = 0.1, rjs = 1.5, r]^ = 7.5, and r]^ = 0.1. The fibers form 
an angle of 7r/4 with the x-axis. 

In Example 2, the initial datum is given by a stimulus applied on the extracellular potential Uq in the 
center of the domain, while both v and the gating variable w are initially set to zero (see Figure 3). The 
units for v, Ue and w are mV. In this example, the following MR setting is chosen. We utilize wavelets with 
r = 3 vanishing moments, a maximal resolution level L = 9, and therefore a finest mesh with Af = 65536 
elements. The reference tolerance given by sr = 5.0 x 10~^. 

We show in Figures 4 and 5 a sequence of snapshots after an initial stimulus applied to the center of the 
domain, corresponding to transmembrane potential v^ extracellular potential Uq and adaptive mesh. 

Table 2 illustrates the efficiency and accuracy of the base MR method by tabulating CPU ratio V, 
compression rate rj and normalized errors. By using MR, we obtain an average data compression rate 
of 17 and an increasing speed-up rate up to 26.09. Moreover, the errors in three different norms remain of 
the order of sr. Here we have computed normalized errors using a reference FV solution on a grid with 
JV = 1024^ = 1048576 control volumes. 

For the time integration using the LTS method, we choose the maximum CFL number allowed by (3.8), 
CFL^=o = 0-5 for the coarsest level and CFL^ = 2^CFL^^o for finer levels. For the RKF computations, we 
use ^desired = 1 X 10""^, Sq = 0.1, iSmin = 0.01, and the initial CFL condition CFLt^o = 0.5. 

We select this example for a detailed comparison of the performance of the FV and MR methods with a 
global time step, the MR method with RKF adaptive global time stepping (MR- RKF), and the MR method 
with local time stepping (MR- LTS). The evolution of the speed-up factor V the and data compression rate 77 
for the MR versions and of the normalized L^ and L^ errors for all these methods are displayed in Figure 7. 
From these plots it is observed that with RKF and LTS, the data compression rate is of the same order 
during the time evolution, which means that the adaptive meshes for both methods should be roughly the 
same. Also, a substantial additional gain is obtained in speed-up rate when comparing with a MR calculation 
using global time stepping: The MR- LTS method gives us an additional speed-up factor of about 2, while 
with the RKF alternative we obtain an additional speed-up of about 4. This effect could be explained in part 
from the lack of need of a synchronization procedure for the RKF computations, and the fact that the CFL 
condition (3.8) is not imposed during the time evolution with the MR- RKF method, allowing larger time 
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Figure 7. Example 2 (bidomain model, one stimulus): Time evolution for data compression 
rate 7^, speed-up rate V, and normalized errors for different methods: MR scheme with global 
time step, MR with locally varying time stepping and MR with RKF time stepping. 



steps. (Although condition (3.8) guarantees numerical stability of the solutions, in practice this is observed 
to be a fairly conservative estimates, and moderately larger time steps may be used.) We can also conclude 
that the errors of the MR-LTS computations are kept of the same order that the errors obtained with a 
global time stepping, while the incurred errors by using the MR-RKF method are larger during the whole 
time evolution. 



6.4. Example 3. For this example, we consider an initial stimulus at the center of the domain, later at 
t = 0.2 ms we apply another instantaneous stimulus to the northwest corner of the domain, and then at 
t = 1.0 ms we apply a third stimulus of the same magnitude to the northeast and southwest corners. The 
system is evolved and we show snapshots of the numerical solution for v, Uq and the adaptive mesh. We use 
the MR-RKF method with JV = 65536, Sr = 2.5 x 10~^, (^desired = 1 x 10~^ and the remaining parameters 
are considered as in Example 2. As in Example 2, from Figures 8 and 9 we clearly notice the anisotropic 
orientation of the fibers. 
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Figure 8. Example 3 (bidomain model, three stimuli): Numerical solution for transmem- 
brane potential v and extracellular potential Uq in [mV], and leaves of the corresponding 
tree data structure at times t = 0.1 ms and t = 0.5 ms. 



7. Conclusions 

We address the application of a MR method for FV schemes combined with LTS and RKF adaptive time 
stepping for solving the bidomain equations. The numerical experiments illustrate that these methods are 
efficient and accurate enough to simulate the electrical activity in myocardial tissue with affordable effort. 
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t = 2.0 ms 

v(x,y,t) 



t = 5.0 ms 

v(x,y,t) 





■ 



I 



I 



Figure 9. Example 3 (bidomain model, three stimuli): Numerical solution for transmem- 
brane potential v and extracellular potential Uq in [mV], and leaves of the corresponding 
tree data structure at times t = 2.0 ms and t = 5.0 ms. 



This is a real advantage in comparison with more involved methods that require large scale computations on 
clusters. We here contribute to the recent work done by several groups in testing whether the combination 
of MR, LTS and RKF strategies is indeed effective for a relevant class of problems. 
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From a numerical point of view, the plateau-like structures, associated with very steep gradients, of typical 
solutions motivate the use of a locally refined adaptive mesh, since we require high resolution near these 
steep gradients only. These areas of strong variation occupy a very reduced part of the entire domain only, 
especially in the case of sharp fronts. Consequently our gain will be less significant in the presence of chaotic 
electrical activity or when multiple waves interact in the considered tissue. 

Based on our numerical examples, we conclude that using a LTS strategy, we obtain a substantial gain in 
CPU time speed-up for a factor of about 2 for larger scales while the errors between the MR- LTS solution and 
a reference solution are of the same order as those of the MR solution. On the other hand, using an MR-RKF 
strategy, we obtain an additional speed-up factor of about 4, but at the price of larger errors. However, in 
assessing our findings, it is important to recognize limitations. The high rates of compression obtained with 
our methods are problem-dependent and they may depend on the proper adjustment of parameters. We 
have only considered here very simple geometries, because all computations are concentrated on adaptivity 
and performance. Simulations on more complex and realistic geometries are part of possible future work. 

Finally, we remark that the FV method given in Section 3 as well the MR framework detailed in Section 4 
are both straightforwardly extensible to the 3D case. A convergence analysis of the implicit version of the 
FV method presented in Section 3 is being prepared [ ]. 
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Appendix 

For illustrative purposes, suppose that we are in the case of a cartesian mesh. We give here an example of 
an interior first-order fiux calculation using LTS for the parabolic part of the bidomain scheme (the equation 
for v), needed to complete a full macro time step, by the following algorithm: 

Algorithm 2 (Locally varying intermediate time stepping). 

(1) Grid adaptation (provided the former sets of leaves). 

(2) do ^ = 1, . . . , 2^ (for the local time steps n + 2"^, n + 2 • 2"^, n + 3 • 2"^, . . . , n + 1^ 

(a) Synchronization: 
do / = L, . . . , 1 

doz = l,...,|AJ,(0. j = l,...,|A|^(/) 
if 1 ^ / ^ Ik-i then 

if V(^ij^^i is a virtual leaf then 

Update reaction terms: 

jn+k2-^ ^ jn+{k-l)2-^ jn+/c2-^ ^ jn+{k-l)2-^ 

'ion{i,j),l ion{i,j),l ' app(i,j),^ app(i,j),^ 

endif 
else 

if V(^j)^^ is a leaf then 

Tn+k2-^ ^ J / n+/c2-^ ^u+k2-^\ 
rn+/c2-^ r / n+fc2-^ n+fc2-^\ 

if V(^+i j)^/ is a leaf then 
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endif 



\^V( 



(i,j),l^^(i,3 + l),l I 



-(^e,(i,j + l),^ - ^e,(i,j),0 



if V(2i+2,2j),^+i; '^(2i+2,2j+i),z+i ^^^ /eai;e5 (interface edges) then 



2(i+l,j),l+l^(2i+l,2j),Z+l 



F^ 



(2i+2,2j + l),^+l^(2i+l,2j+l),^+l 



2(i,j + l),^+l^(2i,2j+l),Z+l 



endif 
endif 
endif 
enddo 
enddo 
(b) Time evolution: 

do / = 1, . . . , L, z = 1, . . . , |AU(0. j = 1, . . . , |A|^(0 
if 1 ^ / ^ Ik-i then t/iere is no evolution: 



F< 



(2i+l,2j+2),^+l^(2i+l,2j+l),^+l 



V 



else 



n+(fc+l)2- 



,,n+/c2- 



Marching formula only for the leaves V(i,j),^.' 

n+(/c+l)2~^ 



V 



_^ n+/c2- 



At; 

Cn 



E 



M^n+/c2-^ 
^ app(i,j),^ 



^^jn+/c2 



(i,j),^ 



I^Vr, 



,^( 



(*,j),i,»/(r.,rr.),Zl 



^(0 <^(^(i,j),^'^(n,m),0 



(^e,(n,m),Z " '^e,{i,j),l) 



endif 
enddo 

(c) Partial grid adaptation each odd intermediate time step: 
do / = L, . . . , /fc + 1 

Projection from the leaves. 
enddo 
do I = liz^ . . . ^ L 

Thresholding, prediction, and addition of the safety zone. 
enddo 
enddo 

Here, Ik denotes the coarsest level containing leaves in the intermediate step /c, h{l) is the mesh size on 
level /, and |A|2(/) is the size of the set formed by leaves and virtual leaves per resolution level / in the 
direction z. The marching formula corresponds to (3.4), for the intermediate time steps /c = 1, . . . ,2^, for 
the leaf V(i,j),^. 

Now we give a brief description of the general multiresolution procedure. 

Algorithm 3 (Multiresolution procedure). 

(1) Initialization of parameters. 

(2) Creation of the initial tree: 

(a) Create the root and compute its cell average value. 

(b) Split the cell, compute the cell average values in the sons and compute details. 

(c) Apply thresholding for the splitting of the new sons. 

(d) Repeat this until all sons have details below the required tolerance £^ 

(3) do n = 1, . . . , total dime .steps 

(a) Determination of the leaves and virtual leaves sets. 

(b) Time evolution with global time step: Compute the discretized space operator A for all the leaves. 

(c) Updating the tree structure: 
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• Recalculate the values on the nodes by projection from the leaves. Compute the details for 
all positions (•, •,/) for I ^ Ik. If the detail in a node and in its brothers is smaller than 
the prescribed tolerance, then the cell and its brothers are deletable. 

• // some node and all its sons are deletable, and the sons are leaves without virtual sons, 
then delete sons. If this node has no sons and it is not deletable and it is not at level 
I = L, then create sons. 

• Update the values in the new sons by prediction from the former leaves. 
enddo 

(4) Output: Save meshes, leaves and cell averages. 

Here total dime steps is the total time steps needed to reach Tfinai using At as the maximum time step 
ahowed by the CFL condition using the finest space step. 

When using a RKF strategy for the time evolution, replace step (3b) by the new step 
(3) • Compute the discretized space operator A for all the leaves as in (5.1). 

• Compute the difference between the two solutions obtained as in (5.2). 

• Apply the limit er for the time step variation and compute the new time step by (5.3). 
When using a LTS strategy, replace step (3b) by the new step 

(3) do n = 1, . . . , total dime _steps 

(a) Determination of the leaves and virtual leaves sets. 

(b) Time evolution with local time stepping: Compute the discretized space operator A for all the 
leaves and virtual leaves 

(c) do A: = 1, . . . , 2^ (k counts intermediate time steps) 

• Compute the intermediate time steps depending on the position of the leaf as explained in 
Section 5.1. 

• if /c is odd then update the tree structure: 

— Recalculate the values on the nodes and the virtual nodes by projection from the 
leaves. Compute the details in the whole tree. If the detail in a node is smaller than 
the prescribed tolerance, then the cell and its brothers are deletable. 

— // some node and all its sons are deletable, and the sons are leaves without virtual 
sons, then delete sons. If this node has no sons and it is not deletable and it is not 
at level I = L, then create sons. 

— Update the values in the new sons by prediction from the former leaves. 
endif 

enddo 

(Now, after 2^ intermediate steps, all the elements are synchronized.) 
enddo 
Here total dim^e steps is the total time steps needed to reach Tfinah with Ato as the maximum time step 
allowed by the CFL condition using the coarsest space step. 
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